* Figure: INTERVENTION STOVE PURCHASE RATES

use "${output}panel_r0_r1_r2.dta", clear // Load merged analysis panel
set seed 293847

********************************************************************************

// Keep households surveyed in each round
keep if sample_household == 1

// Reposition baseline value of household controls
foreach var of varlist household_size children_under_five {
	gen b_`var' = `var' if surveyround == 0
	bys hh_id (surveyround) : replace b_`var' = b_`var'[1] if mi(b_`var')
	}

// Merge with NGO activity raw data
merge m:1 gpname using "${data}chirag_activity_list.dta", ///
	keepusing(chirag_village health education forestry fodder spring youth_club agricultural herbs shg cooperative pirna_ws year_start)

// Sort uniquely for simulations
sort hh_id surveyround

// Keep relevant NGO strata variable
tab chirag_strata chirag_village, miss
drop chirag_village

// Number of simulations
local reps = 10000 // number of bootstrap repetitions

********************************************************************************

* COUNT OF NUMBER OF ACTIVITIES CONDUCTED BY THE NGO

* Generate count NGO-activity intensity count
egen ngo_activity_count = rowtotal(health education forestry fodder spring youth_club agricultural herbs shg cooperative pirna_ws), missing
lab var ngo_activity_count "Count of projects conducted by NGO (maximum possible = 11)"
replace ngo_activity_count = 0 if chirag_strata == 0
gen missing_ngo_activity_count = 0
replace missing_ngo_activity_count = 1 if mi(ngo_activity_count)
lab var missing_ngo_activity_count "=1 if no data available on range of activities"

* Bootstrap with randomly generated NGO activity counts for missing GPs
capture program drop ngo_count_boot_rep
program define ngo_count_boot_rep, rclass
	preserve
		gen ngo_activity_count_replace = runiformint(0, 7)
		replace ngo_activity_count = ngo_activity_count_replace if missing_ngo_activity_count == 1
		bsample, cluster(uniquegrp)
		reghdfe purchased_intervention_stove i.treatment##c.ngo_activity_count b_* if surveyround == 1, absorb(districtcode) vce(cluster uniquegrp)
		return scalar beta_count_interaction = _b[1.treatment#c.ngo_activity_count]
	restore
end

* Do the simulation and generate density plot for full sample
preserve

	qui simulate beta_ngo_count_boot_rep = r(beta_count_interaction), reps(`reps') nolegend seed(249794) : ngo_count_boot_rep

	_pctile beta_ngo_count_boot_rep, percentiles(2.5 5 95 97.5)
	graph twoway (kdensity beta_ngo_count_boot_rep, ///
		xtitle("") ytitle("") ysc(r(0 40)) ylabel(0(10)40)) ///
		|| (kdensity beta_ngo_count_boot_rep, range(`r(r1)' `r(r4)') color(gs14) recast(area)) ///
		|| (kdensity beta_ngo_count_boot_rep, range(`r(r2)' `r(r3)') color(gs11) recast(area) ///
		legend(order(3 2) label(2 "95% C.I.") label(3 "90% C.I.")) name(ngo_count_boot_rep_full, replace) nodraw)  ///
		|| (scatteri 0 0 40 0 , lpattern(dash) lcolor(black) recast(line) title("(a) Count of NGO projects"))

restore

********************************************************************************

* COUNT OF NUMBER OF YEARS OF NGO ACTIVITY

gen years_active = 2013 - year_start
replace years_active = 0 if chirag_strata == 0 
lab var years_active "Number of years NGO has been active in this GP (as of 2013)"
gen missing_years_active = 0
replace missing_years_active = 1 if mi(years_active)
lab var missing_years_active "=1 if no data available on starting year of NGO activities"

* Bootstrap with randomly generated NGO active years for missing GPs
capture program drop ngo_year_boot_rep
program define ngo_year_boot_rep, rclass
	preserve
		gen years_active_replace = runiformint(0, 25)
		replace years_active = years_active_replace if missing_years_active == 1
		bsample, cluster(uniquegrp)
		reghdfe purchased_intervention_stove i.treatment##c.years_active b_* if surveyround == 1, absorb(districtcode) vce(cluster uniquegrp)
		return scalar beta_count_interaction = _b[1.treatment#c.years_active]
	restore
end

* Do the simulation and generate density plot for full sample
preserve	

	qui simulate beta_ngo_year_boot_rep = r(beta_count_interaction), reps(`reps') nolegend seed(24934114) : ngo_year_boot_rep

	_pctile beta_ngo_year_boot_rep, percentiles(2.5 5 95 97.5)
	graph twoway (kdensity beta_ngo_year_boot_rep, ///
		xtitle("") ytitle("")) ///
		|| (kdensity beta_ngo_year_boot_rep, range(`r(r1)' `r(r4)') color(gs14) recast(area)) ///
		|| (kdensity beta_ngo_year_boot_rep, range(`r(r2)' `r(r3)') color(gs11) recast(area) ///
		legend(order(3 2) label(2 "95% CI") label(3 "90% CI")) name(ngo_year_boot_rep_full, replace)) ///
		|| (scatteri 0 0 150 0 , lpattern(dash) lcolor(black) recast(line) title("(b) Years of NGO activity") nodraw)

restore

********************************************************************************

grc1leg ngo_count_boot_rep_full ngo_year_boot_rep_full, ///
	legendfrom(ngo_count_boot_rep_full) cols(1) iscale(*0.9) pos(3) name(allvillages_rep, replace) ///
	b1("{&beta}(TREATMENT{sub:j} × NGO{sub:j})") l1("Density")
	
graph display allvillages_rep, ysize(20) xsize(15)

graph export "${results}figure_ngo_redefinition_linear.png", ///
	width(5000) replace

gr drop _all
